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Abstract 

Pulay's Residual Metric Minimization (RMM) method is one of the 
standard techniques for achieving self consistency in ab initio electronic 
structure calculations. We describe a reformulation of Pulay's RMM 
which guarantees reduction of the residual at each step. The new version 
avoids the use of empirical mixing parameters, and is expected to be more 
robust than the original version. We present practical tests of the new 
method implemented in a standard code based on density-functional the- 
ory (DFT), pseudopotentials, and plane- wave basis sets. The tests show 
improved speed in achieving self consistency for a variety of condensed- 
matter systems. 
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1 Introduction 

The requirement of self consistency between the electronic charge density and 
the potential plays a key role in ab initio electronic structure calculations. The 
iterative search for self consistency generally involves some form of charge den- 
sity mixing at each step. A standard and widely used mixing method is the 
RMM-DIIS technique (Residual Metric Minimization - Direct Inversion of the 
Iterative Subspace), first introducted by Pulay for Hartree-Fock calculations Jl|, 
||, but now also used in density- functional theory ||. The modified Broyden 
technique introduced by Vanderbilt and Louie[Q and SrivastavaQ and gener- 
alised by Johnson[|) can be shown[0, ||] to reduce to the Pulay technique for 
a suitable choice of weights; in practice, these weights give the fastest conver- 
gence, and are often used, though can lead to instability at convergence [|| 
We describe here a new technique for charge-density mixing, and we show that 
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it is more robust and sometimes significantly faster than the Pulay and the 
modified Broyden methods. 

We recall that a self-consistency cycle in an ab initio calculation proceeds as 
follows: an input density p(r) is used to generate a potential v(r); from this, a 
Hamiltonian or Fock matrix is formed; the eigenfunctions of the latter are then 
used to create an output density p'(r). The residual R(r) associated with a given 
input density is defined as R(r) = p'(r) — p(r). The self-consistent density is the 
p(r) for which R(r) = everywhere. The idea of charge mixing in its simplest 
form is that one should attempt to move towards self consistency at each step 
by linearly 'mixing' p(r) and p'(r) in some proportions and using the result as 
input to the next cycle. The Pulay technique is a sophisticated generalization 
of this idea, which works by minimizing the norm | R | of the residual, defined 
as 
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In making a robust method for seeking self consistency, we regard it as highly 
desirable that | R | should decrease at every step. We shall show that, as self 
consistency is approached, our new method achieves this, though the Pulay 
method does not. In the following, we summarize the Pulay method before 
outlining our new formulation, which we call guaranteed-reduction-Pulay (GR- 
Pulay). We then present a number of practical test cases within DFT, which 
demonstrate the advantages of GR-Pulay. 



2 The GR-Pulay method 

At each iterative step, the Pulay method works with a set of s densities, which 
at iteration number n are denoted by p n (r), p„_i(r) . . ., p„_ s+ i(r). (In prac- 
tical calculations, s is typically 5. Of course, in the early stages of the search, 
when n < s, the set of densities is taken to be p„(r),... pi(r).) In going to 
the next iteration, a new density p n+ i(r) is created, and the oldest previous 
density /3„_ s+ i(r) is discarded. The procedure for generating p„ + i(r) involves 
the concept of the present 'optimal' density p° pt (r), which signifies the linear 
combination of the present densities: 

having the smallest norm of the residual, subject to the condition 5Z*_ Q-i = 1. 
To determine the coefficients c^, it is assumed that we are close enough to self 
consistency for variations of the densities and their associated residuals to be 
linearly related. This means that the residual associated with any linear combi- 
nation of the present densities, as in eqn (^J), is simply J2i=o a iRn-i(r), where 
Rn-i(j) is the residual associated with input density p n _^(r). The residual 
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associated with p opt (r), denoted by i?° pt (r), is then given by: 

s-1 

^ pt W = pT(r) ~ P° n pt (r) = £ a z R n ^(r) , (3) 

i=0 

where p n pt '(v) is the output charge density associated with the input p° pt (r). 
Since the residuals R n -i(r) are all known, the constrained minimization of 
J?° pt | is equivalent to the constrained minimization of a bilinear form, and this 
yields a unique and simple formula for the o^. 

The present optimal density p° pt (r) is the 'most self-consistent' linear com- 
bination of present densities. But the new density p n+ i(r) clearly cannot be 
chosen to be p° pt (r). The reason for this is that each new iteration must in- 
troduce new variations of the density, so that p„+i(r) cannot be allowed to lie 
in the subspace spanned by the present densities. In the conventional imple- 
mentation of the Pulay scheme, the new density is usually chosen to be a linear 
combination of p° pt (r) and its output density p° pt '(r): 

Pn+1 (r) = (1 - A)p°?{v) + ApT'iv) ■ (4) 

The value of the mixing coefficient A is empirically chosen (typically to be about 
0.8), and the efficiency of the self-consistency search depends on the choice of A. 
If A is not appropriately chosen, a variety of problems can occur, including slow 
convergence or even failure to converge. The problem is that the best choice of 
A depends on the physical system being treated)^. 

The new mixing scheme proposed here resembles the Pulay method in that 
it works with s densities p n (r), p„_i(r), . . ., p n _ s+ i(r) at each step, and goes 
from one iteration to the next by discarding the oldest density p„_ s+1 (r) and 
creating a new density p„ +1 (r). However, the set of s densities is required to 
have a crucial property: the norm of the residual i? n (r) associated with p„(r) 
is required to be no greater than the norm of the residual associated with any 
linear combination J"^_ a.ip n ^i(r), with the usual condition J"^_ a» = 1. We 
shall see immediately how to ensure this property. Our procedure for generating 
p n+ i(r) is as follows: 

• delete p n - s +i and add p' n , so that the set of densities is p' nl p n , p n -i, ■ ■ ■ 

Pn-s+2- 

• put p' n through a cycle, so that we have its output /?" and hence a full set 
of residuals R' nl R ni R n -i . . . R n - s+2 - 

• make linear combinations: 

p n = a>ip' n + a 2 p n + a3Pn-l + • ■ • + a s p n - s+2 , (5) 

and determine the coefficients ai so as to minimize the norm of the residual 
of p n , subject to the usual condition that the ai sum to unity. Denote by 
p n +i the density p n that yields the minimum residual. 

1 We note that Pulay JlJ |2| implicitly chose A=l; we also note that as the scheme enters the 
linear regime, the value of A becomes unimportant. 
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• delete p' n and replace by p n +i, so that the new set of densities is p n +i, 

Pm ■ ■ ■ Pn-s+2- 

Clearly the new set has the same property as the old: the newest density p n +i 
has the minimal residual norm that can be achieved by taking linear combi- 
nations of the present densities. In particular, | R n +i | < | R n |, so that the 
residual at iteration n + 1 is less than that at iteration n. Moreover, it can be 
shown that the decrease of the residual in an s-level scheme is greater than the 
decrease in an s'-level scheme, provided s > s'. We note that in GR-Pulay, as in 
the original Pulay scheme, the computational effort in each iteration consists of a 
single self-consistent cycle. This is true, provided i? n +i(r) is accurately enough 
given by the linear approximation, which is clearly the case as self-consistency 
is approached. 



3 Applications 

Although our GR-Pulay scheme is completely general, our practical interest here 
is in the pseudopotential, plane-wave implementation of DFT, and our tests 
have been done on a variety of condensed-matter systems using the CASTEP 
codeJlC]. The use of a modified Broyden technique^) (shown by KresseQ and 
Eyertp to reduce to the usual Pulay technique for a suitable choice of the 
weights) in DFT condensed-matter work has been explored by many workers, 
including Johnson and Kresse, and it is one of the standard options in CASTEP; 
it is implemented with the weights chosen so that it is directly equivalent to the 
Pulay technique. The systems chosen for our tests are: bulk silicon (two atom 
primitive cell); bulk magnesium chloride (three atom cell); bulk aluminium (one 
atom stacked fee); bulk plutonium dioxide (three atoms stacked fee); bulk iron 
(one atom stacked fee); and the platinum(OOl) surface (five atoms in five layers, 
slab geometry). These systems provide example of insulating, semiconducting 
and metallic behaviour, as well as extreme inhomogeneity of electron density. 

We have made efforts to achieve an unbiased comparison of Broyden and 
GR-Pulay; all parameters in each run were the same for the two techniques, 
as were the initial conditions. When searching for self-consistency, CASTEP 
applied a criterion of energy change between successive iterations to the Broyden 
method, while we apply a criterion of absolute size of the norm of the residual 
(expressed as a fraction of the norm of the charge density) to the GR-Pulay 
method. We checked that in all cases our criterion was as strict as the energy 
difference criterion (i.e. the norm of the residual when CASTEP had converged 
using Broyden was no smaller than the norm of the residual when CASTEP 
had converged using GR-Pulay). We note that, in the context of the search 
for self-consistency between a charge density and a potential, it is important 
to apply a convergence criterion to the norm of the residual (for this is what 
determines whether or not self-consistency has been reached) and not an energy 
difference; cases where the change in energy is small from one iteration to the 
next, but the norm of the residual is relatively large, can be envisaged. 
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In Table [l] we report timings for the Pulay/Broyden technique and our GR- 
Pulay technique within CASTEP, run on a PC (400MHz Pentium II with 256MB 
of memory) under Linux. The calculations employed standard pseudopotentials 



(either norm-conserving 1 11 |12| or ultrasoft ]T3| ) and plane wave cutoffs of be- 
tween 100 and 300 eV (depending on the system; the efficacy of individual 
methods should be independent of the plane wave cutoff). The criteria for self 
consistency are an energy change of 5x 10 _6 eV per atom for the Pulay/Broyden 
technique and |-R|/|p| < 10~ 4 for the GR-Pulay technique. The first two runs (Si 
and MgC^) involved relaxation of the atomic positions and the third (Al) the 
determination of the unit cell size (these were pursued until the RMS force on 
the atoms was below 0.01 eV/A, and the RMS stress below 0.1 GPa). The next 
three cases (Pu02, Fe and Pt(001)) were all single point energy calculations. 
For the first three, the number of iterations required to reach the self-consistent 
ground state starting from scratch (i.e. at the first iteration) is given in brackets 
after the time. For the last three, the total number of self-consistent iterations is 
given in brackets after the time. Our timings demonstrate that GR-Pulay can be 
over twice as fast as the Pulay/Broyden method; in all cases, this is because GR- 
Pulay requires fewer iterations to find the ground state than Pulay/Broyden. It 
is also related to the number of line searches performed during wave function 
minimisation between self-consistent iterations; modern plane wave DFT codes 
perform several line searches (until either convergence or a maximum number 
of iterations are reached), and differing numbers of these searches will also con- 
tribute to the different times. Monitoring the rate of decrease of the residuals in 
both cases showed that the GR-Pulay method typically achieved a much greater 
rate than Pulay/Broyden (which is not surprising, given the results). 



4 Conclusions 

We have shown that our proposed reformulation of the Pulay mixing scheme 
offers three kinds of advantage over the original scheme. First, no empirically 
chosen mixing parameters are needed, so that less experience is required to 
use the scheme successfully. Second, the residual is guaranteed to decrease at 
every iteration, so that instability cannot occur, at least as self-consistency is 
approached. Third, our practical tests show that the new scheme can sometimes 
be significantly faster than the original one. 
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System 


Pulay/Broydcn 


GR-Pulay 


Si(R) 

MgCl 2 (R) 

Al(CR) 


6.93 (9) 
30.19 (13) 
29.50 (36) 


5.56 (7) 
25.35 (13) 
12.34 (12) 


Pu0 2 

Fe(spin) 

Pt(001) 


20.91 (15) 
63.02 (34) 
670.47 (15) 


16.47 (10) 
33.05 (14) 
688.72 (14) 



Table 1: Timings (in seconds) for the Pulay/Broyden and GR-Pulay methods 
applied to different systems using the CASTEP code. Numbers in brackets after 
the time indicate iterations required to reach the self-consistent ground state 
starting from initial conditions for the first iteration. (R) indicates geometry 
optimisation, (CR) cell optimisation and no letter a single point energy. 
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